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Abstract 

Multicellular tumor spheroids are an important in vitro model of the pre- vascular phase of solid tumors, 
for sizes well below the diagnostic limit: therefore a biophysical model of spheroids has the ability to 
shed light on the internal workings and organization of tumors at a critical phase of their development. 
To this end, we have developed a computer program that integrates the behavior of individual cells and 
their interactions with other cells and the surrounding environment. It is based on a quantitative descrip- 
tion of metabolism, growth, proliferation and death of single tumor cells, and on equations that model 
biochemical and mechanical cell-cell and cell-environment interactions. The program reproduces existing 
experimental data on spheroids, and yields unique views of their microenvironment. Simulations show 
complex internal flows and motions of nutrients, metabolites and cells, that are otherwise unobservable 
with current experimental techniques, and give novel clues on tumor development and strong hints for 
future therapies. 

Author Summary 

When dealing with tumors, a strong emphasis is usually placed on the detailed molecular machinery of 
individual cells. However, their interactions with the environment and their collective behavior are equally 
important, and largely unexplored. Biologists study these interactions in the laboratory, in cultures of 
small spherical clusters of cells called "tumor spheroids" . Unfortunately it is very difficult to exploit the 
full potential of these experimental studies and extract the precious informations withheld in the core of 
spheroids. We have developed a computer model that helps understand the interaction between cells and 
their environment, and establishes a bridge between the microscopic world of molecular interactions and 
the macroscopic properties of spheroids. This computer model is a sort of laboratory where it is possible 
to perform virtual experiments on tumor spheroids and their environment, and gain unhampered access 
to all useful informations on the internal workings of these model tumors. Eventually, we are rewarded 
with novel and unexpected views of the microstructure of tumor spheroids. 

Introduction 

Multicellular tumor spheroids (MTS) stand out as the most important in vitro model of pre- vascular 
solid tumors 1-8 . MTS often have a regular, almost spherical structure, and their apparent simplicity 
has led to repeated attempts to capture their features with neat mathematical models. However, the 
absence of vascularization and the near sphericity hide an internal complexity which is not easy to tame 
either with analytic mathematical models 9-^12^ , or with numerical models based on rough simplifications 
of the biological settings such as cellular automata or other lattice-based models [13 - 16 . Moreover the 
presence of a growing necrotic core 1 1 and of an extracellular matri x pT] , the appearance of convective 
cell motions 18 , and the heterogeneous response to chemotherapies 1 19\ point to the importance of MTS 
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as an in vitro model of tumors, and most of all to their relevance to understand tumor heterogeneity, but 
they also point to the difficulties of producing a useful, predictive model of MTS. 

The appearance of widely different resistance phenomena to antitumor therapies in similarly grown, 
isolated MTS of the same cell type 19 indicates that random fluctuation phenomena play an all-important 
role in the growth kinetics of MTS. It is well-known that the discrete events at the single-cell level (like 
transitions from one cell-cycle phase to the next, mitosis, cell death, etc.) do display some randomness, 
and one can pinpoint the source of large-scale variability on these fluctuations, as they are amplified and 
propagated by cell-cell and cell-environment interactions. Thus, the complexity of MTS development can 
only be captured by a fine-grained, multiscale model, and we need a mathematical description at the 
single-cell level. Since cells communicate with other cells and the environment, the other actors of this 
complex play are the concentration gradients of important molecular species that depend on the structure 
of the extracellular space and of the facilitated transport processes into and out of individual cells, and 
the mechanical forces that push and pull cells as they proliferate with repeated mitoses and then shrink 
after death 20 . These processes mix with complex nonlinear interactions between the biochemical and 
the mechanical part, and this highlights again the importance of an effective model at the single-cell level. 

On the basis of such motivations, we have developed a numerical model of MTS that incorporates 



a working model of single cells 21,22 . We have first put forward a broad outline of its structure in 
reference 23 , and it differs from other models developed in the past [£-16 because it captures at the 
same time both the basic features of cell metabolism, growth, proliferation and death, and provides 
a true lattice- free calculation of cell motions, as they are pushed and pulled by the forces exerted by 
dividing cells, by the growth of other cells, and by the shrinking of dead cells. We also wish to stress 
that the model parameters are either derived from experiment or are deduced from reasonable theoretical 
arguments, so that, essentially, there are no free parameters - there can only be some residual variability in 
biophysically meaningful ranges - the model is truly predictive, and the results are not merely qualitative 
but quantitative as well. 

Here we illustrate in broad terms the structure of the program and report the results of the first 
simulations of single spheroids (technical implementation details are relegated to the supporting text). 
We find that the simulations agree quite well with experimental measurements on real spheroids, and show 
unexpected and important internal patterns. Moreover, we wish to stress that the methods delineated in 
this paper represent very general practical solutions to problems that are common to any simulation of 
cell clusters, and they are just as important. 



Biochemical behavior of individual cells 

The elementary building blocks in this model of MTS are the individual tumor cells that behave as 
partly stochastic automata f2iy22|. Figure [l] summarizes the biochemical pathways that are included in 
the single-cell model: cell metabolism is driven by oxygen, glucose and glutamine, and transforms these 
substances into energy molecules, molecular building blocks and waste products, following the well-known 
biochemical reaction chains f2J. Further details can be found in the original papers 21 22 and in the 



supporting text, which also includes important upgrades to the original model 21,22 . 

In the present version of the program, the stochasticity is mostly concentrated in the discrete events: 
for instance, mitochondria are partitioned at random between daughter cells at mitosis, and cells can die 
because of metabolite accretion, according to a Poissonian cytotoxicity model (see the supporting text). 

We remark that in this approach glutamine also stands for the wider class of aminoacids, and lactate 
is the paradigm of all metabolites: we use the concentrations of glutamine and lactate to represent these 
two classes of substances in phenomenological parameterizations wherever needed. Similarly we use the 
number of mitochondria and ATP content to model the dynamics of cell volume; the single-cell model 
also contains representative members of the cyclin protein class to compute the passage of checkpoints 



and entry into the different cell phases 21 22 25 26 , and finally into mitosis (see also figure SFl in the 
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supporting information for a sketch of the ceh cycle in the simulation program). 

The complete map of the biochemical pathways included in the simulation program is shown in figure 
SF2 in the supporting text. This map comprises only the most basic pathways, however we cannot afford 
to introduce a more complex network at this stage of program development. Indeed, our final aim is the 
simulation of MTS with a volume as large as 1 mm^, which corresponds to more than one million cells, 
so that simulation results overlap actual experimental measurements [19^^27^^ . Since the differential 
system involves 19 independent biochemical variables per cell, we must eventually integrate at least 19 
million coupled nonlinear differential equations for the biochemical cell variables alone (this grows to 
at least 25 million equations when we include the position and velocity variables), and thus even this 
minimal single-cell model leads to a daunting computational task (see the supporting text for further 
details on the algorithmic complexity of the program). 



Reaction-diffusion processes and the environment 

Substances like oxygen are transported into and out of cells by normal diffusion while molecules like 
glucose require facilitated diffusion processes. This means that cell membranes play an important role 
for substances like glucose, and that in this case the diffusion of each such molecular species towards 
cells in the inside of a spheroid needs the free volume in the extracellular space to proceed, and that 
we must model this space as well as the cells to obtain a realistic simulation. We have shown how to 



do this in reference 29 , where we have also discussed ways to tame the exceedingly high stiffness of 
the very large set of reaction-diffusion and transport equations that arise in this context (see also the 
supporting text). The external environment itself is included in these equations, and evolves in time 
as well. In the present model, each cell contributes 15 internal variables and 4 extracellular variables: 
these extracellular variables are the masses of oxygen, glucose, glutamine and lactate in the extracellular 
volume surrounding the cell. Because of its smallness, the extracellular space has an extremely short 
characteristic filling time, which can be as fast as few tens of microseconds. On the other hand, the 
macroscopic features of MTS evolve over times as long as months (i.e., times of the order of lO^s), and 
thus the numerical integrator must be able to handle phenomena that span 12 orders of magnitude in 
time 29 . The internal biochemical reactions included in the numerical model are much slower and their 
fastest characteristic times are only as low as 0.1s, much longer than the diffusion times f29 30 . The 



topology of diffusion in the extracellular spaces is obviously dictated by the cells themselves, and the 
program uses the network of cells centers as the scaffolding for the corresponding discretized diffusion 
problem. The links between the cells' centers - i.e., the proximity relations - are provided by a Delaunay 
triangulation |31l|32| , which is computed repeatedly 33 as the cluster of cells grows and rearranges itself 



under the pushes and pulls of volume growth, mitosis, and the shrinking of dead cells (see also figure 
SF3 in the supporting information). Moreover, the proliferation of cells means that both the number 
of cells and the total number of links steadily grow, and that the differential system of equations that 
model metabolism, transport and diffusion changes all the time, and becomes increasingly complex. The 
3D Delaunay triangulation itself is not an exceedingly heavy computational burden for the program, 
as it turns out that efficient algorithms can compute it, on average, with 0{N) time computational 



complexity 33 35 , so that this algorithm is indeed feasible for very large clusters of cells. 



Biomechanical evolution of the simulated MTS 

Real cells have passive viscoelastic mechanical features, but they also move actively under the pushes 
of their own cytoskeleton, and to the best of our knowledge there is no comprehensive model of cellular 
biomechanics 36,37 . Thus, we resort once again to phenomenological simplifications, and the first and 
foremost is that our cells are stretchable spheres, characterized by their radius, and a few other parameters 
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that specify their viscoelastic properties (see the supporting text for a more detailed description and 
the hst of parameters). We also specify a pairwise interaction force between cells, repulsive when a 
cell pushes against a neighbor, and attractive when we try to detach it from its neighbor. For small 
deviations from the equilibrium distance, we assume that the interaction force is described by the Hertz 
model (explained in the supporting text), while for large deformations due to compression we set the 
force to a fixed saturation value, and for large distances the attractive force decays to zero (see figure 
SF4 in the supporting information) . The description of the interaction forces is tuned to hold also during 
mitosis (see the supporting text and figure SF5). Even though this is a rough approximation of the overall 
mechanical behavior of cells, there are many details that must be managed to make it work, and they 
are all described in the supporting text. 

Here the Delaunay triangulation that we use as the scaffolding for the diffusion problem turns out to 
be useful once again: the same cell-cell links also define the set of neighbors of each cell, and therefore 
the global problem of computing the pairwise interactions between cells can be reduced to a single loop 
over all cells and the small limited number of their immediate neighbors, so that this operation has an 
0{N) computational complexity only - and it does not grow when we include the cost of the Delaunay 
triangulation 35 - instead of the 0(7V^) complexity of generic pairwise interactions. 

Results 

The first and most obvious result is the outstanding match of the growth curves of simulated spheroids 
with those of real spheroids: figure [2] shows a few stages of the growth of a simulated spheroid (a real 
spheroid is shown for comparison in figure [3| , while figure [4] compares the growth curve of a single sim- 
ulated spheroid with the growth curves of real spheroids grown in vitro. Here we see that the growth 
curves are very much alike, and we found that simulation runs with different parameters - in the biophys- 
ically meaningful ranges - produce very similar growth curves, in spite of structural internal changes: the 
growth curves are thus rather robust with respect to parameter changes. 

Several experiments f37||42] have yielded many accurate measurements of oxygen and glucose concen- 
trations and other quantities vs. spheroid radius; these values are part of the output of our simulation 
program as well (see figure [s] and figure [6| , and a comparison with the experimental data is shown in the 
table. On the whole the agreement of simulation data of single spheroids with the experimental values 
is quite good, and we wish to stress that this is not the result of a fit a posteriori, but rather of the a 
priori choice of model features and parameters. These results qualify as true predictions of the numerical 
model. 

The necrotic core of spheroids is another important feature that is well reproduced in the simulations, 
and it is clearly visible in central slices of the simulated spheroid in figure [2] The simulations also provide 
detailed, quantitative snapshots of the necrotic core dynamics; the left column of figure [7| shows the 
percentage of dead cells vs. distance from the centroid of a simulated spheroid at different times. In 
these snapshots we can clearly observe the formation of the sharp step that marks the edge of the necrotic 
core. 

These results indicate that the simulation program is reliable and robust and reproduces - both 
quantitatively and qualitatively - known experimental results. However, it yields much more than just 
successful comparisons: figure [8] shows two views of the spheroid microenvironment that at present would 
be unobtainable by other means at this level of resolution. The left panel of figure [S] is a plot of the fiow 
of glucose in the extracellular spaces of a mature spheroid, superposed on a density plot of extracellular 
glucose concentration, and it shows - rather unexpectedly - that there is an outward fiow of extracellular 
glucose from the central necrotic region. In the external, viable rim the fiow is inward bound, and there 
is a spherical shell where the flow is stationary. The right panel of figure [8] shows the corresponding 
plot of cell velocities in the same central slice, and we see that the velocity vectors point outward in the 
viable rim, while there are well-formed vortices in the central region, and the region in-between displays 
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distinctive chaotic motions: these three regions closely match the three regions in the left panel. The 
right column in figure [7| shows radial velocity vs. distance from the centroid of the simulated spheroid, 
and sheds some more light on the nature of this structure: as more and more cells die and the necrotic 
core forms, the dead cells shrink and the core contracts. The contraction of the necrotic core expels the 
residual glucose in the extracellular spaces and produces the observed outward fiow. We found that this 
behavior is strongly dependent on the particular value of the diffusion coefficient and on the metabolic 
activity of cells. In some simulations - where we used a lower value for the effective diffusion coefficient 
of oxygen - we observed a similar structure with oxygen as well. We remark that in the case of lactate we 
found no such structure, and we obtained a pH value - derived from the distribution of lactate inside the 
spheroid - that is very close to experimental measurements: this indicates that the discretized reaction- 
diffusion scheme used in the simulation program performs correctly, and that the observed flows are not 
algorithmic artifacts. 

Discussion 

Although the program described in this paper is based on a model of individual cells that includes 
only the basic cell functions, the simulation results compare very well with experimental measurements, 
and give strong hints on the sources of individual spheroid variability. Moreover, the images obtained 
in single runs reveal unexpected and interesting correlations and an elaborate structure of the tumor 
microenvironment that could never be observed before. This unexpected, complex microstructure - the 
formation of different regions, and the ffows that characterize them, along with the complex velocity 
field - can be discerned in the flows of the other substances, though not all of them, according to their 
effective diffusion coefficient and their metabolism: the figures of these fiows are shown at full-size as 
supporting information. Thus if we suppose that, in a more complete description, there are N substances 
that characterize the spheroid microenvironment, and assume that the spherical shell that divides the 
two main regions lies in the same position for all of these substances and that their effective diffusion 
coefficients are uncorrelated, then 2^ different spheroid structures are determined by diffusion alone. 
The variation of some critical parameter (e.g., a slight change in the metabolic activity due to local 
fluctuations in the number of dead cells, and thus a change in the effective diffusion coefficients) can 
potentially act as a switch and determine widely different fates for similar spheroids. This variability 
cannot be discerned from growth experiments: the simulations that we have performed to date indicate 
that the growth curve alone is not enough to distinguish between such different states, because it does 
not change much even when important substances, like oxygen, diffuse in markedly different ways. These 
different states represent different biochemical conflgurations of tumor microenvironment, that might 
exert distinct selective pressures on cells during tumor evolution. 

The spheroid microstructure that is well evidenced in flgurejsj and in flgures SF8-21 and in the movie 
flies Sl-3 in the supporting information, shows highly correlated fluctuations that produce, e.g., islets 
of proliferating cells in the sea of dead cells of the core, and cell and mass flows that follow preferential 
channels. There is a sort of spheroid-speciflc self-organization of the internal structure due to these 
correlated fluctuations. Similar cell flows have been observed in the lab and a recent review has stressed 
the great signiflcance of such flndings 43 : the simulations suggest that the whole topic of cell flows 
and extracellular diffusion should be investigated further. On the basis of the simulation results, we also 
conjecture that the ffow of therapeutic drugs may be diverted as well, and let some viable, proliferating 
tumor cells escape treatment. This means that the simulation program could eventually become an 
important tool to design novel treatment schedules, and possibly validate the effects of anti-tumor drugs. 

Certainly the model is far from complete, and we plan to add soon several new features, like a basic 
model of intracellular acidity, now accounted for by a simple phenomenological parameterization, and the 
effects of pH and salt concentration on diffusion. However, already in its present form, we believe that 
this numerical model is a true testbed of biological complexity and a real virtual laboratory, and also a 
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source of important biomedical clues. 



Methods 

The simulation program is written in ANSI C++: this programming language was a natural choice from 
the very start for distinct reasons: 

• C++ is an object-oriented language, and in a simulation such as this, it is very natural to define 
objects that have a clear-cut biological meaning; 

• at present, C++ programming is supported by a vast array of scientific libraries, and this helps 
reducing program development time; 

• the availability of the fiexible and powerful C++ library CGAL 33 that handles the computational 
geometry structures utilized by the program (convex hulls, Delaunay triangulations and alpha 
shapes); 

• the availability of powerful development tools and highly optimized compilers. 

The structure of the program refiects the organization explained in the paper: a layout is shown in 
figure |9] The functional blocks work as follows: 



Initialization 

At start, the internal variables of all cells are set at approximate standard values. During initialization, 
cells are allowed to grow and proliferate freely in an environment that is held fixed. The number of 
cells is also kept constant, and when a mitosis occurs one of the daughter cells is discarded. In this 
initial phase cells can have large oscillations of their metabolic parameters, and can occasionally step in 
parameter regions that would normally spell death: this does not occur here. Initialization lasts until 
the oscillations of metabolic parameters die out. We have determined the duration of the initialization 
phase observing the desynchronization of a population of initially synchronized cells: when oscillations 
of the relative fractions of cells in each cell-cycle phase become undetectable we estimate that cells have 
reached a stable state. It turns out that a simulated time of 3 • lO^s (i.e. about 35 days of simulated 
time) is sufficient for initialization of cell with a period of about 20 hours. Usually the starting number 
of cells is quite small (normally just one cell to seed the growth of a single spheroid), and initialization 
executes in very short real time (a few seconds). 



Metabolism, diffusion, transport, and growth 

This part of the program solves the combined differential system of equations that describe internal cell 
metabolism and diffusion in the extracellular spaces (described in detail in the supporting text), using 
the implicit Euler method. This leads to a system of nonlinear equations, that are solved in turn with a 
variant of the Newton-Raphson method. The functional scheme of this important part of the program is 



shown in figure 10 We wish to stress that although the number of variables can be quite large (more than 
10^ loop variables), convergence is reasonably fast, because the initial concentration values are invariably 
very close to the final ones. 



Cell motion 

Cell motion is also regulated by differential equations and the solution uses a strategy based on a semi- 
implicit method (described in detail in the supporting text). Volume growth is regulated by the part 
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that handles metabohsm and diffusion, therefore it is loosely coupled to cell motion. However we have 
implemented an updating mechanism that effectively decouples the two parts of the program: this means 
that the program can use multithreading with shared memory and exploit the features of multicore 
processors. 



Cellular events 

This part of the program handles discrete events, like cell-cycle transitions, mitosis and cell death. In 
case of mitosis it also initializes the daughter cells - using the metabolic variables of the mother cell - 
and allocates memory for the new cells. 



Geometry and topology of cell cluster 

Geometrical and topological informations are updated here, using calls to CGAL methods 33 that 
compute the convex hull of the cluster of cells, the Delaunay triangulation of cell centers, and the alpha 
shape of the cluster - with an alpha parameter [33 equal to (2ro)^ where ro is the average cell radius. This 
part of the program uses this basic information to set all relevant geometrical and topological variables 
in the program. 



Summary statistics and dump on file 

The last step in the loop computes several statistics and outputs them on summary files. It also writes 
periodically the whole configuration of cells on file for further processing. 



Program termination 

The program repeats the loop until one of the stop conditions is met: either all cells are dead, or the 
program executed the required number of steps. The supporting information text contains additional 
considerations on algorithmic complexity and on measured performance (see figure SF6 and figure SF7). 

Additional processing to extract useful informations from the simulation data is performed with several 
standard tools, like Mathematica f44] 



Supporting Information 

Supporting text: provides technical details on the structure of the simulation program and includes 
tables STl to ST5. 

Figure SFl: sketch of the cell phases included in the simulation program. 
Figure SF2: sketch of the metabolic network. 
Figure SF3: the geometry and topology of diffusion. 

Figure SF4: pictorial representation of the interaction force between two cells. 
Figure SF5: the geometry of mitosis. 

Figure SF6: CPU time needed to simulate 1 hour, vs. the number of cells in the spheroid. 

Figure SF7 : total CPU time vs. the total number of cells. 

Figures SF8-10: oxygen concentration. 

Figures SFl 1-14: extracellular glucose concentration. 

Figures SF15-17: extracellular glutamine concentration. 

Figure SF18: lactate concentration . 

Figures SF19-21: velocity in the plane of the slice. 

Movie SI: development of the necrotic core. 
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Movie S2: flow of extracellular glucose. 
Movie S3: map of projected cell velocities. 
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Figures 




Figure 1. Rough sketch of the biochemical pathways incorporated in the model of single 
cells. We take into account the main metabolic pathways (glycolysis, oxidative phosphorylation 
through the TCA cycle and gluconeogenesis), including the role of mitochondria in the production of 
ATP. The model also includes protein and DNA synthesis, and checkpoints controlled by representative 
members of the cyclin family The single-cell model has two spatial compartments (the inside of the cell 
and its immediate neighborhood, the extracellular space that surrounds it) and transport of substances 
between these compartments is regulated by transporters on the cell membrane that are also included 
in the model. The extracellular space of each cell communicates by simple diffusion with the 
neighboring extracellular spaces and with the environment. The complete map of the biochemical 
pathways is shown in figure SF2 in the supporting information text. 
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Figure 2. Snapshots of one simulated spheroid taken at different times. As the spheroid 
grows, a necrotic core develops in its central region, just as it happens in real spheroids. The size of the 
necrotic core and of the viable cell rim match real measurements. 




100 |jm 



Figure 3. Photograph of a spheroid grown in vitro from HeLa cells in agar. The spheroid is 
colored with trypan blue to mark dead cells, where the necrotic core is clearly visible. The agar 
contains the spheroid and helps in obtaining a better spherical shape with HeLa cells, but also stifles 
spheroid growth because it reduces the effective diffusion coefficients in the nourishing medium, so that 
it cannot be directly compared to the simulated spheroid in the second column of figure [2] (which has 
the same size), while it is similar to the larger spheroid in third column. 
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Figure 4. Growth curve of a simulated tumor spheroid (solid line). The run parameters used 
in this case are hsted in tables ST3, ST4 and ST5 in the supporting information text. The symbols 
denote data points taken in different in vitro experiments: squares = FSA cells 
(methylcholantrene-transformed mouse fibroblasts) 45 ; diamonds = MCF7 cells (human breast 
carcinoma) [19^; circles = 9L cells (rat glioblastoma) |27j 
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Figure 5. Concentrations in the simulated spheroid. The color coded figures on the left show 
the partial pressure of oxygen, the concentrations of glucose and lactate in the extracellular spaces, and 
the pH of the extracellular environment (high values = red, low values = blue). The corresponding 
plots in the right column show the average values of these quantities vs. the distance from the centroid 
of the tumor spheroid. The small oscillations in the plots close to the spheroid surface are due to 
fluctuations in the averaging procedure, because the spheroid is slightly nonspherical. 
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Figure 6. Plots of the normalized average intracellular concentration of lactate (green), 
glucose (blue), and ATP (red). These plots have been obtained in the same simulation and at the 
same time step as the plots of figure |5) and each concentration is normalized to its peak value. These 
plots indicate that cell death in the central region is due both to the accumulation of metabolites 
(lactate) and to metabolic stress (starvation). 
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Figure 7. Fraction of dead cells (left column) and average radial velocity (right column) at 
different times. As the spheroid grows, the necrotic core becomes increasingly weh defined, and as 
dead ceUs shrink, the radial velocity changes sign and a marked inward motion characterizes the central 
region. 
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Figure 8. Two views of the microstructure of a simulated spheroid, with about 500/im 
diameter and 296264 cells (183893 live cells + 112371 dead cells). (Left panel): flow of 
extracellular glucose along a central section of the tumor spheroid (yellow arrows) superposed on the 
plot of glucose concentration. The length of the arrows is proportional to the glucose flow intensity 
projected on the plane of the section. At this stage, the necrotic core is contracting as dead cells 
gradually shrink, and this leads to a slow outward flow of the glucose stored in the extracellular spaces 
in this central region. We observe that such a behavior depends on the effective diffusion coefficient of 
glucose, and it disappears completely when the diffusion coefficient is high enough. This also suggests 
that the flow of glucose and other substances, like therapeutic drugs, is strongly dependent on the 
biochemistry and structure of extracellular spaces, and even small changes can lead to markedly 
different internal spheroid morphologies. (Right panel): individual cell velocities in the simulated 
spheroid. This is the same central section as in the left panel, and the velocity vectors are projected on 
the plane of the section. The length of each vector is proportional to the projected speed. The velocities 
in the viable rim show a coherent outward motion, while the velocities in the necrotic core show a 
rather orderly inward motion, with some vortices due to local residual cell proliferation. The region 
in-between is somewhat chaotic and the global structure of this plot mirrors that of the glucose flow 
shown on the left. The supporting information includes higher-quality versions of these figures and 
those of other flows. 



Model of Tumor Spheroids 



18 



Initialization 



loop 

(fixed timestep) 



Metabolism, diffusion, 
transport, and growth 



Cell motion 



Cellular events (cell cycle, 
mitosis, and death) 



Geometry and topology of cell 
cluster 



Summary statistics and dump 
on file 



Figure 9. Functional blocks of the simulation program. Program initialization is followed by a 
loop that performs biochemical and biomechanical calculations. This is followed by a check of the status 
of individual cells - this is where we decide whether a cell advances in the cell cycle, undergoes mitosis, 
or dies. Next the program computes the geometry and the topology of the cell cluster, and finally it 
outputs intermediate statistics and results. The loop continues until a user-defined stop condition is 
met. Some parts of the program can proceed in parallel (like metabolism and cell motion), and we can 
use multithreaded code. 
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Figure 10. Functional blocks of the C++ method that computes metabolic and 
extracellular variables. This part performs a loop that computes the solution of the nonlinear 
equations found in the implicit Euler integration step [29] (see also the supporting text). Although the 
number of variables can be quite large (more than 10^ variables), convergence is fast, because the initial 
concentration values are invariably very close to the final ones. 
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Table 1. Comparisons with experimental parameters 



Parameter 


Simulation 


Experiments 


^Cell type 


Ref. 


^Glucose uptake (kg m~'^) 


1.44-10-^ 


5.4-12.6-10-^ 


Rat-Tl, MRl 


f37 


"'^ Lactate release (kg m~'^) 


1.35 • 10"^ 


5.4 - 9 • 10"^ 


Rat-Tl, MRl 


[37 


^p02 (mmHg) 


7 


0-20 


Rat-Tl 


[37 






0-40 


MRl 


37 






20-60 


EMT6 /Ro 


38, 




6.7 


6.6 


C6, H35 


|39 






6.96-6.99 


U118-MG, HTh7 




40. 


^ApH 


0.77 


0.41 


U118-MG 




4o| 






0.49 dz 0.08 


HTh7 




40 


^Viable cell rim thickness (/im) 


155 


200 


EMT6/R0 




38 






142 ± 16 


HTh7 




40 






310 ± 28 


U118-MG 




40 






198 ± 27 


C0II2 




41 








225 ± 26 


HT29 




41 




^Hypoxic rim thickness (/im) 


98 


44 ± 52 


C0II2 




41 








44 ± 52 


HT29 




41 




Cell cycle distribution (%) 


GO/Gl = 57.3 


GO/Gl = 58 ± 4 


BMG-1 




42. 




S = 21.6 


S = 19 ± 1 








G2/M = 21.1 


G2/M = 23 ± 1 







Metabolic and histologic parameters in spheroids of approximately 500 /im diameter: comparison 
between a single, large simulation, carried out with the parameters listed in tables ST4 and ST5 in the 
supporting text, and experimental data. 

Notes: 

1. Rate of glucose uptake or lactate release per viable spheroid volume (see f37]). 

2. Central p02 tension (experiments) or estimated in the centroid (simulations). 

3. pH has been determined in the central region of the spheroids. This corresponds to a sphere of 
radius ~ 100/im about the centroid of the spheroid. 

4. Difference between environmental pH and pH 200 jam below the spheroid surface. 

5. In our simulations the viable cell rim thickness corresponds to the distance between the spheroid 
surface and the inner shell where only 5% of the cells are still alive. Experimental values have been 
determined by histology. 

6. These values corresponds to the radius of the necrotic core. 

7. Cell types are as follows: Rat-Tl = T24Ha-ras-transfected Rati cells (Rati = spontaneously im- 
mortalized rat embryo fibroblasts); MRl = myc/ T24Ha-ras-cotransfected rat embryo fibroblasts; 
EMT6/R0 = mouse mammary tumor cells; C6 = rat glioma cells; H35 = rat hepatoma cells; 
U118-MG = human glioblastoma cells; HTh7 = human tyroid carcinoma cells; Coll2 = moderately 
differentiated human colon adenocarcinoma ; HT29 = poorly differentiated human colon adenocar- 
cinoma; BMG-1 = human glioma cells. 



